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ABSTRACT 


“=. Introduction 

More than two-thirds of the area of Lebanon are prone to karstification processes including surface karst features 
such as poljes, uvalas, dolines, blind valleys, natural bridges, as well as smaller features like karren [1]. 
Subsurface karst features comprise many types of solutional shafts and galleries, caves, subsurface lakes and rivers 
and most types of speleothems. The Jurassic, Cretaceous and Eocene limestones /dolomites form the core of the 
Lebanon ranges [2]. Sinkholes play a major role in recharging amounts of groundwater, mainly from the melting 
snow covering the Lower Jurassic and Upper Cretaceous (Cenomanian) rocks several months a year. The 
groundwater recharge process in the area benefits as well from rainfall events [1]. Sinkholes and closed 
depressions are used as detention basins and ponds since decades and, thus, are important for the water supply and 


water management. 


The development of karst phenomena, the intensity and depth of karstification are influenced by various 
geodynamic movement such as uplift/subsidence or folding, by the structural pattern and by volcanic activities [3]. 
The pattern of zones of secondary porosity such as faults and fracture zones are controlling factors of the 


groundwater flow and, thus, of the karst development in the investigation area. 
1.1. Study Objectives 


This study aims to contribute to the mapping of the surface karst features, especially sinkholes, and to identify their 


orientation, concentrations and dimensions, as well as the tectonic pattern and further factors influencing their 
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development, based on remote sensing and GIS derived data and references. A systematic GIS approach aims to 


gather those factors, to study their complex interactions, and to create a GIS data base. 


Of course, the evaluation capabilities of satellite data are limited, and additional research is required to clarify the 
character and properties of the mapped karst features. This is a requirement for land use planning as the collapse of 
dolines or caves might be a risk for settlements and the infrastructure. Therefore, one goal of this study is to 


provide hints and point out directions where future field research or other surveys could be focused. 


: 2. Materials and Methods 


This study is dealing with the evaluation of optical satellite data as well as of radar images and of Digital Elevation 
Model (DEM) data and DEM derived morphometric maps, aiming to detect surface traces of karst features, fault 


zones and their related structures, combined with other geoscientific data of Central-Lebanon. The evaluations of 


different satellite data and available geodetic and geophysical data were carried out as summarized in Figure 1. 
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Figure 1. Materials, software, and workflow. 


Various factors and parameters derived from thematic layers gained from satellite data and references such as 
lithology, structural lineaments, terrain analysis and earthquakes were combined in a GIS data base to achieve 


additional knowledge or to contribute to a better understanding of the known facts. 
2.1. Optical Satellite Data 


Digital image processing of Landsat TM and 8/9 (the Operational Land Imager (OLD, Sentinel 2 and ASTER data 


was carried out by merging different Red Green Blue (RGB) band combinations. Special attention was focused on 
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the creation of RGB images including Landsat 8/9 thermal bands 7 and 10 and ASTER thermal bands (10-14). 
Filter tools were used for the structural analysis of Landsat and Sentinel 2 data. Sentinel 2 data were included into 
the data sets. As digital image processing software served the Sentinel Application Platform (SNAP)/ESA and 
ENVI/L3Harris Geospatial Solutions as well as the processing tools integrated into the geoinformation systems 
ArcGIS/ESRI and QGIS. The specific digital image processing methods of the different satellite data and image 
enhancements methods are focused on gaining additional knowledge about the structural pattern such as the filter 
techniques. The data were provided by the USGS Earth Explorer [4], the ESA Copernicus Browser [5], and the 
NASA Earth Data [6]. 


The high-resolution satellite images such as OrbView-3 images, Bing Map from Microsoft and Google Earth 
images as well the World Imagery data from ArcGIS Earth from ESRI were used as base for detailed mapping. 


2.2. Radar Data 


Sentinel 1 Synthetic Aperture Radar (SAR, C-Band radar wavelength), ALOS PALSAR (L-Band, 23.5 cm 
wavelength) radar data were digitally processed and evaluated. The Advanced Land Observing Satellite (ALOS), 
Phased Array type L-band Synthetic Aperture Radar (PALSAR) from the Japan Aerospace Exploration Agency 
(JAXA) [7] provided data with 12.5 m spatial resolution. The Sentinel 1 and ALOS PALSAR radar missions 
included dual polarization capabilities. The radar system can transmit a signal in a horizontal (h) or vertical (v) 
polarization, and then receive in both h and v polarizations. RGB- tools were used to merge radar images with 
different polarizations and, thus, get “false-colored” radar images improving the evaluation possibilities, for 
example by Sentinel 1-vh, vv, vh polarization composites. The different radar wavelengths correlated with the 


L-Band (23.5 cm) and C-Band (5.6 cm) allow varying radar penetration depths into the surface. 


The illumination geometry of the radar signals plays an important role for the visibility of the linear and structural 
features, especially in the Central-Lebanon area with its high mountain ranges oriented in SSW-NNE direction, 
creating effects as well like radar image distortions like layover and radar foreshortening. Those oriented 
perpendicular to the radar illumination angle appear enhanced, whereas those parallel to the radar illumination are 
more suppressed, return large amounts of the transmitted energy to the radar satellite sensors (visible in light gray 
tones), while smooth surfaces scatter the energy and, thus, have low radar returns such as the Bekaa valley bottom 
(appearing dark on the radar images). Other factors such as dielectric properties and radar characteristics 
(frequency, depression angle, polarization) also affect the radar return [8], [9]. The processing of the radar data was 


carried out using the SNAP software of ESA and the image processing tools integrated in ArcGIS. 


The radar data were provided by the ESA Copernicus Browser [5], and NASA Earth Data, Alaska Satellite Facility 
(ASF) [10], [11]. ALOS PALSAR mosaics (25 m-resolution) were retrieved from the Earth Observation Research 
Center, PALSAR global mosaic and Forest/Non-forest Map, JAXA [7]. 


2.3. DEM Data 


Digital Elevation Model (DEM) data (SRTM, ASTER and ALOS PALSAR) and the DEM derived morphometric 


maps were included into these investigations in a Geolnformation System (GIS) embedded environment. 
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Morphometric maps like height level and slope gradient maps support the detection of deformation structures 
along major fault zones. The position of slopes oriented towards wind-driven, rain-bearing clouds has an impact on 
the intensity of karst development. Therefore, aspect maps derived from DEM data were used to provide 


information about the orientation of slopes. 


For the detection of areas with a relatively higher surface water input the Hydrology tools of ESRI were used to 
calculate the drainage system (stream order) and flow accumulation as well as the drainage basins based on the 
DEM data. Weighted overlay tools were used to find areas where the susceptibility to higher surface water input 
(up to flash floods) is higher due to the accumulation of specific morphometric properties such as the lowest local 


height level, low slope gradients, no curvature of the terrain or high flow accumulation values. 
2.4. Structural Analysis 


Shapefiles were created in ArcGIS related to karst features, surface properties, faults and to structural features. 
Based on the different satellite data a visual lineament and structural analysis was carried out. The visibility of 
linear features depends on the specific properties of the satellite systems such as their spatial and radiometric 
resolution, but also on the acquisition time. The seasonal and climatic conditions, as for example the phenological 
situation of the vegetation, soil moisture conditions, sun elevation, illumination angle and direction at the 
acquisition time of the satellite data play an important role. The detailed inventory of the surface-near fault and 
fracture pattern helps to detect those areas, where the permeability of the lime- and dolostones is relatively higher, 
supporting intrusion and infiltration of meteorological surface water after heavy rains. Special attention was 
directed at distinct expressed linear features (tonal linear anomalies, geomorphologic linear features, drainage 
segments, etc.). In the scope of this study the following types of linear and curvi-linear features were mapped: 
lineaments (as a neutral term for linear features without precisely knowing their origin, fault zones, and structural 
features. The visual lineament analysis was chosen because it allows the differentiation of several types of linear 
features what is not possible when using automatic lineament extraction in this high mountain environment. The 
extracted linear geological features were imported into a GIS environment and statistical analysis and 
representation was performed using rose diagrams. Density calculations were performed. The linear traces of 
faults were merged in a GIS with available geological and geophysical information (geological and geophysical 
maps georeferenced into ArcGIS and QGIS) such as the active faults from the Global Active Faults Database [12]. 


Infrastructural data were downloaded from Geofabrik [13]. 


Earthquakes can trigger the development of fault and fracture zones, caves, dolines and large conduits collapse, 
and change groundwater flows, not only because of the earthquake shock, but also because of horizontal and 
vertical movements. To learn about areas prone to relatively higher ground motion, earthquake data were 


downloaded from different sources [14], [15], [16], and the epicenters integrated into the GIS data base. 


:3. Geographic and Geologic Overview 


As karstification processes in Central-Lebanon are strongly connected to the geologic and geographic conditions a 


short overview of the environmental setting is given. 
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3.1. Geographic Setting 


The Central-Lebanon is characterized by two mountain ranges that run parallel to the coast of the Mediterranean 
Sea. The Mount Lebanon range borders a narrow coastal plain. The Mount Lebanon and Anti-Lebanon ranges are 
separated by the Bekaa valley (Figure 2a and Figure 2b). The Bekaa valley formed by a graben structure, is an 
elongated valley located between the Mt. Lebanon and Anti-Lebanon ranges. The NNE-SSW trending mountain 
ranges with their steep slopes on both sides and the Beqaa valley between the ranges, as well as the narrow costal 
plane are all imprints of the structural deformation during the geologic past. This complex relief results in a high 
spatial variability of precipitation [17] which is increasing due to climate change [18]. High intensity, but short 
rainfall events cause flooding, even flash floods. High precipitation is documented just to the north of Beirut with 


1600 to 1700 mm of precipitation per year! [19]. 
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Figure 2a. Height level map based on SRTM DEM data. 
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Figure 2b. Topographic cross section in W-E-direction of the Central-Lebanon. 
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3.2. Geologic Overview 


The present structural framework of Lebanon is mainly shaped by the activity of the Dead Sea Transform fault 
(DSTF) that connects the ongoing opening of the Red Sea and Gulf of Aqaba to the Tauros-Zagros subduction 
zone. A dominant element is the Bekaa Graben with the left-lateral Yammouneh shear fault at its western 
boundary [2]. The Bekaa Graben is considered to be the northern extension of the Red Sea-Dead Sea-Jordan Rift 
Valley transform fault, the Dead Sea Transform Fault. The structural setup is mainly governed by elements 
resulting from the African/Arabian-European continental collision between the Upper Cretaceous and the Tertiary. 
Along the Yammouneh fault average late Pleistocene slip rates of 5.1 + 1.3 mm/year were monitored based on the 


offset of 25-ka-old alluvial fans [20]. 


The geological units occurring in the study area are mainly composed of limestones /dolomites, volcanic rocks and 
clastic sediments. Gravity data suggest that the Beqaa valley is filled by up to 9 km of lacustrine and continental 
sediments with the deepest part in the central area of the valley. The timescale of the limestone formations ranges 
from the Lower Jurassic to Upper Cretaceous. Unroofed Jurassic units along the entire length of the Mt. Lebanon 
range suggests laterally constant uplift [20]. Bends and stepovers and pull-apart depressions along the strike-slip 
fault segments related to the Yammouneh fault trace stress heterogeneities, affecting the regional strain field. 
Geometrical relationships between folds and strike-slip features suggest that regional strain partitioning may 


accommodate a convergent component of motion between the Arabian and African plates [21]. 
= 4, Results 


According to the specific properties of the remote sensing data the gained results were subdivided. Each satellite 
system offers different advantages (and disadvantages) regarding the evaluation results. Because of the 


combination of the evaluation results additional knowledge can be gained. 
4.1. Evaluations of optical satellite data 


High resolution satellite images (Bing Map, World Imagery, OrbView, Google Earth) served as base for the 
digitizing of surface karst features like sinkholes, enclosed depressions and poljes. Sinkholes and smaller enclosed 
topographic depressions are often used as detention and retention ponds. This is verified by comparative 
evaluations of older high resolution satellite images such as provided by Google Earth (Historical Imagery tool) 
showing sinkholes and closed, bowl-shaped depressions, that were later modified to detention ponds for water 
management. However, the mapping of sinkholes can lead to errors whenever detention basins were created 
without a sinkhole existence before. (A detention basin does not have a permanent pool of water. The water runs 


out or evaporates between rainfalls and, thus, it often remains dry). 


Nevertheless, all the water-filled ponds and circular depressions visible on the satellite images were mapped as 
karst related depressions, assuming that the circular detention basins of the same size as sinkholes were built using 
the existing, natural depressions (Figure 3). In the coast-near, densely populated, sealed areas many sinkholes are 
not visible anymore on satellite images. Thus, errors related to the sinkhole-detection and mapping based on 


remote sensing data cannot be excluded. 
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Based on the sinkhole-shapefiles density calculation were carried out and information of the sinkhole’s size 
derived. Most of the digitized sinkholes have small sizes between 10-60 m (Figure 4). For the better visibility 


buffer zones of 300 m were drawn around the sinkholes. 
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Figure 3. Sinkholes modified to detention basins to preserve water availability as visible on Bing Map images. 
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Figure 4. Length/diameter of sinkholes and detention basins and ponds. 


When analyzing the distribution pattern of the sinkholes in Central-Lebanon, a nearly semi-circular pattern 
becomes visible, concentrated in a belt-like, nearly semi-circular zone, surrounding the highly fractured 
Cretaceous (Cenoman) limestones with brittle surfaces. This structure is bordered in the eastern part by the large 
Yammouneh fault zone. Most of the visible sinkholes are situated in the areas between 1500 and 2000 m height 
level. The areas above 2000 m height level are covered by snow for several months per year, thus, providing melt 
water during the warmer months. Infiltration of snow melt and surface water run-off intensifies the karstification in 


those areas. 


Traces of major fault zones, their parallel segments, and fault related structures like pull-apart-depressions forming 
poljes (-flat-floored depression within karst limestones), smaller push-up ridges, and brittle surfaces are clearly 
visible on the different satellite images that could be verified as such by the comparison with available geologic 
and tectonic maps (Figure 5 and Figure 6) [22], [23], [2], [3], [24]. Understanding the impact of deformation 


patterns on the hydraulic properties of the rock is essential for a better knowledge of karst processes. Uplifting 
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movements, intrusions of volcanic rocks, earthquake shock, and high tectonic stress have contributed to the 


development of the brittle surfaces supporting dissolution processes and collapses in the lime- and dolostones. 


Brittle surfaces are characterized by the rough surfaces including small, parallel valleys and ridges, shallow, 
bowl-shaped depressions and sinkholes, often in a linear arrangement (uvalas). Karstification processes contribute 
to the sculpturing and modelling of the brittle surfaces by creating the small enclosed depressions, poljes, 


sinkholes and karren surfaces. 


The brittle surfaces are concentrated along the SSW-NNE striking fault zones. The brittle surfaces can be observed 
along the SSW-NNE oriented major fault zones, the Yammouneh strike-slip fault, of the Bekaa valley and their 
intersections with SW-NE striking faults (Figure 6). An overview of the digitized fault zones, brittle surfaces and 
karst features is shown in Figure 5 and Figure 6. Various types of deformation due to geodynamic activity could be 


observed (Figure 7). 


Sinkhole 
Buffer_Sinkholes_300 m 
a 

AFEAD_v2022-Faults 
Major fault zone 
Lineaments 


Structural features 
Brittle surfaces 
Pull-apart depression 
a 

Volcanic rocks 
gis_osm_waterways 
gis_osm_water 


MEDKAM _Karst_v1_poly 


MEDKAMclas, MEDKAMclas, MEDKAMclas 
[1 Karst aquifers in evaporite rocks, Karst aquifers in evaporite rocks 
[Bl Karst aquifers in sedimentary and metamorphic carbonate rocks 
|" Various hydrogeological settings in other strata 


Figure 5. 3D perspective view of the structural pattern of Central-Lebanon looking towards north. 


The information of the volcanic rocks war derived from the geologic map [22]. Known faults (brown lines) were 


included from [12] and [25]. 


As there exist outcrops of volcanic intrusions in this area (Figure 5), it has still to be further investigated to which 
extent the uprise of magmatic bodies underneath and intrusion into the limestones and dolomites has caused 
thermal altering, metamorphism and fracturing. Magmatic intrusions (here mainly basalts) can form a hindrance 


and natural barriers for the groundwater flow in the affected areas, especially when occurring as dykes. 


The degree of brittleness in the limestones is mainly controlled by their mineralogic composition. It increases with 
higher quartz and dolomite content, whereas an increase in the clay content represents an increase in ductility [26]. 
The brittle-ductile transition depends on grain size and porosity. Porous carbonate rocks are often mechanically 
weak, heterogeneous at multiple scales, sensitive to effective stress changes, and to the chemistry of the fluid 
saturating the pore space. (-Brittle deformation happens when rocks fail as rigid solids. The rocks will break, rather 
than bend, under these conditions to produce fractures. Brittle deformation occurs along discrete planes in the rock 
instead of involving the rock body as a whole. Later, fractures can be places where material is removed as a result 


of rock dissolution, or minerals can grow in the open fracture spaces.) [27]. The several months long duration of 
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snow cover and, thus, later the melting water contributes to dissolution processes and to the modeling the brittle 
surfaces. Due to the dense drainage pattern, surface water is “trapped” in the small, narrow valleys. They function 
obviously to funnel rainfall and downwards runoff and, thus creating rillen-karren. The dissolution conditions in 
the limestones within the small channels support the development of closed, bowl-shaped depressions, often 


arranged as uvalas, and of elongated poljes. 
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Figure 6. Structural evaluation of satellite data and overview of brittle surfaces in Central-Lebanon. 
The green line indicates the nearly semi-circular outline of the distribution of sinkholes. 


Faults from the Global Active Faults Database [12] were included. 
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Figure 7. Types of deformation monitored by the evaluation of different satellite data. 
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The brittle surfaces vary in their brittle density as shown in Figure 8 (areas 1-3). The drainage pattern allows the 
detection of the structural setting. The orientation of the surface fault and fracture pattern is tracing the stress 


situation changes. 


The small, narrow valleys and ridges are clearly visible on the Modified Normalized Difference Water Index 
(MNDWIJ) image result created with the SNAP software (Figure 9a). The MNDWI image reveals information 
related to soil moisture, in this case of snow cover. The small valleys and depressions, still covered by snow, 


appear in blue colors. The red lines represent the Yammouneh shear fault segments. 
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Figure 8. Brittle surfaces of outcropping limestones with different brittle intensity patterns 1- lower density, 


2 and 3 — higher density) intersected by faults (arrow). 


The drainage pattern and orientation reflect the stress situation clearly [28]. The majority of the small segments are 
oriented parallel to the SSW-NNE striking Yammouneh shear fault. The rose diagram of Figure 9c, derived from 
the drainage orientation analysis based on ALOS PALSAR DEM data (calculated in QGIS), shows the main 


directions of the drainage segment orientations. 


The dense fracture pattern facilitates the infiltration of surface water and, thus, karstification processes such as 
towards karren karst forms (rill-like erosional form of limestone solution). Sinkholes and springs appear to be 
concentrated at the western border of the nearly semi-circled outline of the limestones with brittle surfaces in the 
north of the city of Zahle (Figure 10). The information of the springs was provided by [3] and [25]. Their 
distribution is influenced by SW-NE and SSW-NNE striking fault zones. The distribution pattern of sinkholes is 
also an important factor for the water recharge, especially for those sinkholes forming uvalas along predominant 
alignments. These alignments are usually tracing faults that create high permeability and porosity zones in the 


affected areas [1]. 
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Figure 9a. Brittle surfaces visible on a Modified Normalized Difference Water Index (MNDWI). 
MNDWI= (Green - SWIR)/(Green + SWIR) image, color-coded, derived from the 
Sentinel 2 satellite data acquired on 24.04.2024. 


Figure 9. b) Bing Map scene of a brittle surface and c) drainage pattern tracing the structural pattern 


and karren landforms and the derived drainage segment orientation (trend shown in red colors) 


in the west of the Yammouneh fault zone (red lines). 
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Figure 10. Distribution of springs according to [3] and [25]. 


Sinkholes occur concentrated within a nearly semi-circular belt-zone. 
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Earthquake epicenters were included into the investigation [14,15,16] and their density calculated (Figure 11 a and 
b). When analyzing the earthquake pattern during the two last decades it becomes obvious that stronger ground 
motions after earthquakes must have contributed to the intense geomechanically stress within the limestones 
resulting in a dense fracture and deformation pattern and, obviously, as well as to long term horizontal and vertical 


movements. Most of the earthquakes happen in depths below 10 km. 
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Figure 11b. Earthquake density calculation. (a) Distribution of sinkholes, brittle surfaces and earthquakes. 
(b) Amplification of (a). Sinkholes (red) occur concentrated along the intersection of 


SW-NE and SSW-NNE striking fault zones with a higher earthquake density. 


The thermal bands of Landsat satellites support the detection of circular features (Figure 12). Their origin might be 
related to larger sinkholes or to local uprise of magmatic bodies up-doming the strata above. Outcrops of 
Cretaceous and Tertiary volcanic rocks (mainly basalts, often as dyke intrusion) are documented on the geologic 
maps of this region, occurring in the “belt zone” of sinkhole concentration. In case of volcanic influence, rising 
pressure and temperature due to limestone—magma interaction might be the reason for the denser brittleness. 
Hydrothermal activity related to volcanism with increased dissolution obviously led to the concentration of karst 


features near volcanic intrusions. Traces of ring structures as visible on the Landsat scene (Figure 13, including in 
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the RGB image the thermal bands 7 and 10) might be explained by circular and ellipsoidal zones around collapse 
structures after dissolution processes in limestones as well. The corresponding Google Earth scene of the circular 
structure is presented in Figure 14. Water-filled, small circular depressions are visible as dark spots. Enclosed, 
bowl-shaped depressions show partly a linear arrangement. Not all of them seem to be sinkholes. To clarify their 


properties field investigations are necessary. 
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Figure 12. Structural features and sinkhole distribution concentrated in a nearly semi-circlar pattern visible on a 
Landsat 8 scene of the central Bekaa valley. The position of the circular structure shown in 


Figure13 and Figure 14 is indicated by the black rectangle. 


35S130E 35°52'0°E 35°52'30°E 
1 05 0 1 Kilometers 
es 


RGB2740-LC08_L1TP_174036_20231105_20231110_02_T1 


Figure 13. Circular structure visible on a Landsat 8-scene (30 m spatial resolution) situated within 


limestones with brittle surfaces. 
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Figure 14. 3D perspective Google Earth view (position: 34°02'22.97" N 35°52'10.15" E) of the 


central part of the circular structure, partly covered by snow. 
4.2. Evaluations of Radar Data 


Although radar shadow, foreshortening and lay-over effects on the radar images have to be taken into account 
because of the high mountain terrain conditions, some additional information can be derived as demonstrated in 
Figure 15 and Figure 16. Closed circular depressions can be detected on the L-Band images of the ALOS PALSAR 


satellite that do not appear on optical satellite images (Figure 16, white arrows). 


35°S10°E 38°520°E 35°530°E 35°54'0°E 35°5S'0'E 35°56.0E 35°57 0'E 
2 1 ° 2 Kilometers 


RGBHV,HH,HV-ALPSRP243860670-H2.2_UA. 


Figure 15. ALOS PALSAR radar scene tracing circular features. 


HV-ALPSRP243860670-H2.2 UA 


Figure 16. Radar scene (ALOS PALSAR) revealing smaller sinkholes (black spots indicated by white arrows) and 


fracture zones (dark lines, indicated by red arrows). 
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4.3. Evaluations of DEM Data 


SRTM, ASTER (both about 30 m spatial resolution) and ALOS PALSAR (about 12.5 m resolution) DEM Data 
were processed in ArcGIS using the Spatial Analyst tools to derive different morphometric maps, among them the 
delineation of drainage basins and morphologic watersheds (Figure 17). However, as karst catchment areas 
comprise a complex water transport system characterized by the heterogeneity of surface and underground karst 


forms, this drainage basin map is covering only the information of surface water-runoff. 


The drainage pattern (stream order according to Strahler) was derived as well from the DEM data (Figure 18). To 
get an idea where the highest surface water input can be expected, the drainage density was calculated based on the 
drainage shapefiles derived from ALOS PALSAR DEM data (Figure 18). (- The analysis of drainage network 
density represents the sum of fluvial segment length inside a grid cell). Higher drainage densities can be observed 
along the larger fault zones, especially along the Yammouneh shear fault, and along a nearly semi-circular zone 
with a width of about 5 — 10 km in the west of the Bekaa valley, as indicated in Figure 18b. The visible sinkholes 
occur concentrated within the same zone. The highest concentrations of the drainage segments are visible on 
Figures 18 a and b in the northeast of the city of Zahle. The higher the meteorological surface water input, the more 


intense get the karstification processes and the higher is the groundwater recharge potential. 


BSISWE IS4OVE ISKSTE 'SOUE ISSUE 36°00 ESO IOUE 36ISOE 36°200E 


70 5 0 10 Kilometers 
as 


Figure 18a. Drainage density calculation derived from ALOS PALSAR DEM data. 
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Figure 18b. Drainage density calculation merged with the occurrence of visible sinkholes. 


Some morphometric factors influence the disposition to be prone to relatively more intense karstification such as 
height levels, slope orientations, slope degree, and curvature of the terrain, etc. A weighted overlay procedure can 
be carried out for the detection of areas with higher susceptibility to flash floods and to higher surface water input 
and infiltration by extracting first causal/preparatory factors and, then, by aggregating these morphometric factors 
(lowest local height level, slope degrees < 5°, curvature = 0, aspect = (-1) flat) in the weighted overlay-tool of 
ArcGIS (Figure 19). The result of the aggregation of these morphometric factors shows areas susceptible to higher 
surface water input after precipitations and even susceptible to flash floods due to their morphologic disposition. 
This approach is useful as well to detect pull-apart depressions, poljes and the morphometric properties of the 


larger valleys. 
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Figure 19. Weighted overlay approach in ArcGIS - Susceptibility to higher surface water input and infiltration (4 


-7) due to the local morphometric disposition. 


The exposition of the slopes plays an important role when considering the main rain bearing air masses, clouds and 
wind directions, here coming from SW and W. West winds arriving from the Mediterranean Sea are often changed 
in their directions towards north or south due to the hindrance of the high mountain ranges, the Mount Lebanon 


ranges and Anti Lebanon ranges, forming a wind barrier (Figure 20). Because of the orographic situation and 


ISSN: 2581-5059 OPEN @ACCESS 148 


v 


oe ne: Volume 8, Issue 2, Pages 133-155, April-June 2024 


S Mediterranean Journal of Basic and Applied Sciences (MJBAS) 


effects the slopes oriented towards SW get more rain (Figure 21) as air masses are forced to flow over the high 
topography. Whenever air rises over the Mount Lebanon ranges, it cools, and water vapor condenses. As a result, 
the rain occurs concentrated on the windward side of the mountain range. Rainfall increases with elevation in the 
direction of wind directions and cloud tracks. Sinkholes in the limestones occur exactly where the highest 


precipitations were measured (Figure 22). 


The majority of the sinkholes are situated luv-wards in the west of the highest areas of the Mount Lebanon range. 
Most of the visible sinkholes are situated along slopes oriented towards SW to W what can be verified by aspect 
maps. By calculating aspect maps from DEM data, the slopes oriented toward W and SW can be derived. Thus, the 
relationship between the exposition of slopes to rain-bearing clouds and relatively higher surface water input 


(leading to more intense karstification) and the occurrence of sinkholes can be visualized. 


Orographic Effects on Rain bearing Clouds 


Wind Streaming Pattern https://earth.nullschool.net 


Figure 20. Wind pattern (19.03.2024) indicating the orographic effects of the Lebanon mountains by directing the 


western winds towards north. Wind data were downloaded from [29]. 


Orographic effects on rainfall distribution in Central-Lebanon Fill-SRTM-Mosaic 


m 

Ho- 100 
Increasing rainfall 
with elevation 


3D perspective view based on SRTM DEM data 


180° 
Wind frequency rose diagram— Mount Lebanon 


Figure 21. Orographic effects on the precipitation amount (schematically represented by blue arrows) 


considering the main wind directions. 
Rose diagram: https://globalwindatlas.info/en/area/Lebanon/Mount%20Lebanon 
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Figure 22. Annual precipitations according to WorldClim version 2.1 climate data for 1970-2000 
released in 2020 [19] with 30 seconds (~1 km2) spatial resolution. 


The orientation of slopes derived from DEM data by calculating aspect maps is shown in Figure 23. 
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Figure 23b. Aspect-map combined with the mapped sinkholes and a typical cloud situation image 
(acquired: 22.11.2022) [30]. 
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Pull-apart depressions and basins become visible when subtracting the original DEM data from the Fill-DEM data 
created in ArcGIS using the tools in the Spatial Analyst extension. The first step used the “Fillsink” algorithm from 
the ArcMap software (identifying the point or set of adjacent points surrounded by neighbors with higher elevation 
and rises to the lowest value on the depressions boundary). This procedure then fills all depressions in the DEM, 
including both those generated from data errors (spurious artifacts) and those that record real topographic features, 
such as karst depressions (dolines). The second step was to extract the sink depths in these areas by differencing 
the maps between the sink-filled (“depressionless” DEM) and original DEM. The difference image can highlight 
the different depressions, including the karst enclosed depressions. Figure 24 shows the result of such a subtraction 
of SRTM-DEM from Fill-SRTM DEM data. The elongated pull-apart-depressions along the left-lateral shear zone 
appear very precisely. This method can be applied to detect larger sinkholes and poljes as well whenever the sizes 
depths of the bowl-shaped and elongated depressions are distinct enough to be visible (Figure 25). However, the 
DEM spatial resolution (SRTM and ASTER DEM - 30 m, ALOS PALSAR DEM - 12.5 m) is not sufficient 
enough to detect the smaller sinkholes and enclosed depressions. Smaller depression cannot be detected as 
demonstrated by the Google Earth scene with depressions that are not visible on the calculated difference map 
(Figure 25). 


Figure 24. Detection of pull-apart depressions (blue) based on ALOS PALSAR DEM data 
using the difference between (Fill-DEM — DEM). 


Figure 25. Detection of larger depressions (blue) by calculating (ALOS PALSAR_Fill DEM — DEM). 
Smaller depressions along the narrow valleys and ridges are not visible 


due to the limited spatial resolution of the DEM data. 
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5. Conclusion 


The multisource approach presented in this article contributes to a further knowledge of the development of karst 
in this region. The different satellite and geodata support the better understanding of the occurrence and 
distribution pattern of surface karst features in Central-Lebanon as far as they are visible on the satellite images. 
The use of morphometric factors, especially the spatial distribution of the drainage network, drainage density 
calculations, orientation analysis, combined with tectonic lineaments data and lithologic layers provide additional 
knowledge related to the hydrogeological characterization of the region. The evaluations of the different satellite 
data reveal the existence of structural features such as ring structures, that, sometimes, become only visible on the 


thermal bands of the Landsat images or on radar images. 


Several factors seem to have an influence on the spatial distribution and concentration of karst features as 
summarized in Figure 26: Occurrence of limestones, heights > 1500 m, slopes oriented towards W, SW, NW, 
slopes with slope degrees between 40-70°, annual precipitations above 1000 mm/year, larger fault and fracture 
zones (especially when intersecting each other) and a higher earthquake density. The secondary effects of 
earthquakes in this area and their characteristics have to be considered, depending on their magnitudes, depths and 
spatial concentrations. The predominantly shallow earthquakes below 10 km depth with higher magnitudes have 
an impact on the rock properties, especially on the development of fractures and faults. Thus, in areas with a higher 
earthquake density more surface karst features can be observed. Larger fault and fracture zones influence the 
hydrogeologic conditions due to the higher water permeability in the affected areas. The analyses of the different 
satellite images and of DEM data contribute to the detection of areas with higher surface water input which is 


important for the groundwater recharge. 


) The occurrence of the majority of surface karst features in Central- 
Lebanon is influenced by the following environmental, geomorpho- 
logical, lithological and geophysical settings: 

* Limestones, dolostones 


* Heights > 1500m 
* Slopes oriented towards W, SW, NW 


* Slopes with slope degrees between 40-70° 


+ Annual precipitations between 1000— 
1600 mm / year 


+ Fault and fracture zones 


* Earthquake density, magnitude and depth 


Figure 26. Main factors influencing the occurrence of sinkholes (brown spots) 


and other surface karst features in Central-Lebanon. 


As karst water resources are important for the water supply and economic development in Lebanon, more detailed 


research of the specific local situation in the karst areas is carried on. The structural evaluations of remote sensing 
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data can support the detection of fault and fracture zones in karst areas with influence on groundwater 
permeability. However, high resolution image and DEM data are required for a more detailed analysis such as 
gained by drones or by the LIDAR techniques. They are necessary as well for the detection of sinkholes. Future 
research should include the investigation of the impact of magmatic intrusions into the limestones in Central 
Lebanon and, as consequence, on aquifers. Long term horizontal and vertical movements play an important role, 
too, affecting the development of fracture zones. More geodetic data have to be considered. Further on, the effects 
of climate changes with a relatively higher frequency of extreme weather events have to be monitored in the karst 


areas such as extreme draughts or flash floods and their impact on the karst aquifers. 
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